
clc;
clear;

% Load point estimates (Y) and estimate of variance matrix (Sigma)
regression_output = readtable('RegressionOutput.csv');

regression_output = table2array(regression_output);

Y = -regression_output(:,1:3)';

Sigma = [regression_output(:,4), regression_output(:,7), regression_output(:,8);...
         regression_output(:,7), regression_output(:,5), regression_output(:,9);...
         regression_output(:,8), regression_output(:,9), regression_output(:,6)];

% Set table = 3 for Table 3 and table ~= 3 for Table 8
table = 3;

% 'Nominal significance level'
alpha = 0.05;
     
if table == 3
    % In Table 3, we assume that the researcher is willing to impose a priori
    % knowledge that the two treatment effects (of "therapy only" and "cash
    % only") are non-negative 
    bounded = [1;1;0];
else
    % In Table 8, we assume that the researcher is willing to impose a priori
    % knowledge that all three treatment effects are non-negative 
    bounded = [1;1;1];
end

bounds = [0;0;0];

% 'standard errors'
std_errors = sqrt(diag(Sigma));
     
tool = 1./sqrt(diag(Sigma));
     
% 'Correlation matrix'
Omega = diag(tool) * Sigma * diag(tool);

%% Upper one-sided confidence interval for "therapy only" 
'Row 1'

SSCI_TO = SSCI(Y,Sigma,1,1,alpha,bounded,bounds);

SCI_TO = [Y(1)-norminv(1-alpha)*std_errors(1),Inf];

excess_length_SSCI_TO = Y(1)-SSCI_TO(1);

excess_length_SCI_TO = Y(1)-SCI_TO(1);

ratio_TO = excess_length_SSCI_TO/excess_length_SCI_TO;

fprintf('$')
fprintf('%.4f', Y(1)) 
fprintf('$ & $')
fprintf('%.4f', std_errors(1)) 
fprintf('$ & $[')
fprintf('%.4f', SSCI_TO(1)) 
fprintf(',\\infty)$ & $')
fprintf('%.4f', excess_length_SSCI_TO)
fprintf('$ & $[')
fprintf('%.4f', SCI_TO(1)) 
fprintf(',\\infty)$ & $')
fprintf('%.4f', excess_length_SCI_TO)
fprintf('$ & $')
fprintf('%.4f', ratio_TO)
fprintf('$\n')

%% Upper one-sided confidence interval for "cash only"
'Row 2'

SSCI_CO = SSCI(Y,Sigma,2,1,alpha,bounded,bounds);

SCI_CO = [Y(2)-norminv(1-alpha)*std_errors(2),Inf];

excess_length_SSCI_CO = Y(2)-SSCI_CO(1);

excess_length_SCI_CO = Y(2)-SCI_CO(1);

ratio_CO = excess_length_SSCI_CO/excess_length_SCI_CO;

fprintf('$')
fprintf('%.4f', Y(2)) 
fprintf('$ & $')
fprintf('%.4f', std_errors(2)) 
fprintf('$ & $[')
fprintf('%.4f', SSCI_CO(1)) 
fprintf(',\\infty)$ & $')
fprintf('%.4f', excess_length_SSCI_CO)
fprintf('$ & $[')
fprintf('%.4f', SCI_CO(1)) 
fprintf(',\\infty)$ & $')
fprintf('%.4f', excess_length_SCI_CO)
fprintf('$ & $')
fprintf('%.4f', ratio_CO)
fprintf('$\n')

if table == 3
    %% Two-sided confidence interval for "both" 
    'Row 3'

    SSCI_B = SSCI(Y,Sigma,3,0,alpha,bounded,bounds);

    SCI_B = [Y(3)-norminv(1-alpha/2)*std_errors(3),Y(3)+norminv(1-alpha/2)*std_errors(3)];

    length_SSCI_B = SSCI_B(2)-SSCI_B(1);

    length_SCI_B = SCI_B(2)-SCI_B(1);

    ratio_B = length_SSCI_B/length_SCI_B;

    fprintf('$')
    fprintf('%.4f', Y(3)) 
    fprintf('$ & $')
    fprintf('%.4f', std_errors(3)) 
    fprintf('$ & $[')
    fprintf('%.4f', SSCI_B(1)) 
    fprintf(',')
    fprintf('%.4f', SSCI_B(2))
    fprintf(']$ & $')
    fprintf('%.4f', length_SSCI_B)
    fprintf('$ & $[')
    fprintf('%.4f', SCI_B(1)) 
    fprintf(',')
    fprintf('%.4f', SCI_B(2)) 
    fprintf(']$ & $')
    fprintf('%.4f', length_SCI_B)
    fprintf('$ & $')
    fprintf('%.4f', ratio_B)
    fprintf('$\n')

    %% Two-sided confidence interval for "interaction effect" 
    'Row 4'

    R = [1 0 0;0 1 0;-1 -1 1];

    % Note that now point estimates for "therapy only", "cash only", and 
    % "interaction effect" are stored in Y
    Y = R*Y;

    Sigma = R*Sigma*R';
    
    % 'standard errors'
    std_errors = sqrt(diag(Sigma));

    tool = 1./sqrt(diag(Sigma));

    % 'Correlation matrix'
    Omega = diag(tool) * Sigma * diag(tool);
    
    SSCI_I = SSCI(Y,Sigma,3,0,alpha,bounded,bounds);

    SCI_I = [Y(3)-norminv(1-alpha/2)*std_errors(3),Y(3)+norminv(1-alpha/2)*std_errors(3)];

    length_SSCI_I = SSCI_I(2)-SSCI_I(1);

    length_SCI_I = SCI_I(2)-SCI_I(1);

    ratio_I = length_SSCI_I/length_SCI_I;

    fprintf('$')
    fprintf('%.4f', Y(3)) 
    fprintf('$ & $')
    fprintf('%.4f', std_errors(3)) 
    fprintf('$ & $[')
    fprintf('%.4f', SSCI_I(1)) 
    fprintf(',')
    fprintf('%.4f', SSCI_I(2))
    fprintf(']$ & $')
    fprintf('%.4f', length_SSCI_I)
    fprintf('$ & $[')
    fprintf('%.4f', SCI_I(1)) 
    fprintf(',')
    fprintf('%.4f', SCI_I(2)) 
    fprintf(']$ & $')
    fprintf('%.4f', length_SCI_I)
    fprintf('$ & $')
    fprintf('%.4f', ratio_I)
    fprintf('$\n')
else
    %% Upper one-sided confidence interval for "both" 
    'Row 3'

    SSCI_B = SSCI(Y,Sigma,3,1,alpha,bounded,bounds);

    SCI_B = [Y(3)-norminv(1-alpha)*std_errors(3),Inf];

    excess_length_SSCI_B = Y(3)-SSCI_B(1);

    excess_length_SCI_B = Y(3)-SCI_B(1);

    ratio_B = excess_length_SSCI_B/excess_length_SCI_B;

    fprintf('$')
    fprintf('%.4f', Y(3)) 
    fprintf('$ & $')
    fprintf('%.4f', std_errors(3)) 
    fprintf('$ & $[')
    fprintf('%.4f', SSCI_B(1)) 
    fprintf(',\\infty)$ & $')
    fprintf('%.4f', excess_length_SSCI_B)
    fprintf('$ & $[')
    fprintf('%.4f', SCI_B(1)) 
    fprintf(',\\infty)$ & $')
    fprintf('%.4f', excess_length_SCI_B)
    fprintf('$ & $')
    fprintf('%.4f', ratio_B)
    fprintf('$\n')
end


